clear 
set more off

cap log close
log using "${logs}/theoretical_bias.log", replace

local h_ADGP1 = 0.334449
local h_BDGP2 = 0.173211
local h_CDGP3 = 0.238465
local h_DDGP4 = 0.31935
local h_EDGP5 = 0.147428
local h_FDGP6 = 0.288458
local h_GDGP7 = 0.159048
local h_HDGP8 = 0.214066
local h_IDGP9 = 0.101346
local h_JDGP10 = 0.148273
local h_KDGP11 = 0.190064

local tau_lq_ADGP1 = 9.809037
local tau_rq_ADGP1 = 6.924716

local tau_lq_BDGP2 = 96.22992
local tau_rq_BDGP2 = -113.4154

local tau_lq_CDGP3 = -3.556184
local tau_rq_CDGP3 = -6.357568

local tau_lq_DDGP4 = -1.288178
local tau_rq_DDGP4 = -3.678939

local tau_lq_EDGP5 = -15.81248
local tau_rq_EDGP5 = 15.99251

local tau_lq_FDGP6 = 0.4000707
local tau_rq_FDGP6 = -0.7679749

local tau_lq_GDGP7 = -8.654692
local tau_rq_GDGP7 = 1.088098

local tau_lq_HDGP8 = -1.789442
local tau_rq_HDGP8 = 0.1467769

local tau_lq_IDGP9 = -21.28845
local tau_rq_IDGP9 = 18.03547

local tau_lq_JDGP10 = 3.12233
local tau_rq_JDGP10 = -0.3569933

local tau_lq_KDGP11 = 2.586736
local tau_rq_KDGP11 = -0.2259673


foreach paper in "ADGP1" "BDGP2" "CDGP3" "DDGP4" "EDGP5" "FDGP6" ///
			"GDGP7" "HDGP8" "IDGP9" "JDGP10" "KDGP11" {	

	use "${dat}/DGP/`paper'_ScaledRunningVariable.dta", clear
		
	local c = 0
	local p = 1
	local deriv = 0
	
	local h = `h_`paper''
	
	** tau_rq and tau_lq are the second-order derivatives, as opposed to the regression coefficient
	** on the second-order term
	local tau_lq = `tau_lq_`paper'' * 2
	local tau_rq = `tau_rq_`paper'' * 2
	
	tempvar x x_l x_r
	
	qui gen `x' = x
	qui gen `x_l' = `x' if `x'<`c'
	qui gen `x_r' = `x' if `x'>=`c'	

	qui su `x_l'
	local N_l = r(N)

	qui su `x_r' 
	local N_r = r(N)

	local N = `N_r' + `N_l'

	
	mata {
	
	X   = st_data(.,("`x'"),   0)
	X_l = st_data(.,("`x_l'"), 0)
	X_r = st_data(.,("`x_r'"), 0)
	c = `c'
	N_l = length(X_l)
	N_r = length(X_r)
	c_l  = J(N_l, 1, c)
	c_r  = J(N_r, 1, c)
	p1 = `p' + 1	

	wh_l = kweight(X_l,`c',`h',"triangular");	wh_r = kweight(X_r,`c',`h',"triangular")
	uh_l = (X_l-c_l)/`h';	uh_r = (X_r-c_r)/`h';
	uhh_l=select(uh_l,wh_l:> 0);	uhh_r=select(uh_r,wh_r:> 0)
	
	Xh_l  = select(X_l,  wh_l:> 0);	Xh_r  = select(X_r,  wh_r:> 0)
	whh_l = select(wh_l, wh_l:> 0);	whh_r = select(wh_r, wh_r:> 0)
	Nh_l = length(Xh_l); Nh_r = length(Xh_r)
	X_lp  = J(N_l,p1,.);	X_rp  = J(N_r,p1,.)	
	Xh_lp = J(Nh_l,p1,.);	Xh_rp = J(Nh_r,p1,.)
	
	for (j=1; j<=p1; j++)  {
		X_lp[.,j]  = (X_l:-c):^(j-1);		X_rp[.,j]  = (X_r:-c):^(j-1)
		Xh_lp[.,j] = (Xh_l:-c):^(j-1);		Xh_rp[.,j] = (Xh_r:-c):^(j-1)
							}	
		
	Hp_vec = J(p1, 1, .)
	for (j=1; j<=p1; j++) {
		Hp_vec[j] = `h'^(-(j-1))
							}
	Hp = diag(Hp_vec)	
	
	Gamma_lp  = quadcross(X_lp,wh_l,X_lp);	Gamma_rp  = quadcross(X_rp,wh_r,X_rp)	
	invGamma_lp  = cholinv(Gamma_lp);	invGamma_rp  = cholinv(Gamma_rp)	
	
	v_lp = (Xh_lp:*whh_l)'*(uhh_l:^(`p'+1));	v_rp = (Xh_rp:*whh_r)'*(uhh_r:^(`p'+1)) 
	
	BiasConst_lp = factorial(`deriv')*cholinv(Hp)*invGamma_lp*v_lp 
	BiasConst_rp = factorial(`deriv')*cholinv(Hp)*invGamma_rp*v_rp
	Bias_tau = (`tau_rq'*BiasConst_rp[`deriv'+1,1] - `tau_lq'*BiasConst_lp[`deriv'+1,1])*(`h'^(`p'+1-`deriv')/factorial(`p'+1))
	
	st_numscalar("Bias_`paper'", Bias_tau)
			}
			
	di "Asymptotic bias of local linear estimator for `paper': " Bias_`paper'		
				}


log close
